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ABSTRACT 

The  complex  behavior  and  failure  mechanisms  of  fiber-reinforced  composite 
materials  have  made  multiscale  modeling  a  necessity  when  determining  the  response 
of  composite  structures.  Current  multiscale  models,  which  bridge  information  between 
the  microstructural  constituents  and  the  macroscopic  response  of  a  composite  material, 
are  limited  to  a  prescribed  set  of  boundary  conditions  to  kinematically  link  the  macro 
and  micro  scales,  with  the  most  popular  being  periodic.  This  restriction  prevents  these 
types  of  models  from  accounting  for  the  effects  of  locally  non-periodic  regions  within 
a  composite  structure.  One  such  region  is  a  free-edge  boundary,  which  is  a  common 
damage  initiation  zone  for  many  composite  structures.  In  this  work,  a  semi-concurrent 
multiscale  model  is  implemented  within  ABAQUS  which  allows  for  non-periodic 
boundary  conditions  through  the  development  of  an  energy  based  constitutive 
coupling  between  the  scales  whereby  Hill’s  energy  condition,  or  energetic  consistency 
between  the  macro  and  micro  scales,  is  preserved.  The  methodology  is  applied  to 
examine  free-edge  effects  on  2D  lamina  RVE’s  with  cohesive  fiber/matrix  failure 
within  a  semi-concurrent  scheme. 

INTRODUCTION 

The  complex  and  hierarchical  nature  of  composite  materials  make  them  highly 
desirable  for  high-performance  and  multi-functional  applications.  Consequently,  the 
various  length  scales  and  reinforcement  architectures  present  in  a  given  composite 
pose  a  serious  challenge  for  engineers  and  designers  in  predicting  and  modeling 
failure  and  damaged  response  [1].  Over  the  last  several  decades,  a  variety  of 
multiscale  methods  have  been  developed  to  incorporate  the  effects  of  the  composite 
microstructure  on  the  overall  macroscopic  behavior  [2],  Available  multiscale  models 
can  be  classified  into  three  generic  categories:  sequential,  concurrent,  and  semi¬ 
concurrent  modeling  [3].  These  classifications  are  presented  schematically  in  Figure  1. 

Sequential  models  include  standard  homogenization  schemes  which  determine 
macroscopic  material  properties  based  on  the  composite  microstructure  and 
microconstituent  properties  [4-7],  These  models,  however,  do  not  preserve 
microstructural  information  post-homogenization  and  only  provide  initial  elastic 
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Composite  Vehicle  Research  Center,  Michigan  State  University,  2727  Alliance  Dr.,  Lansing, 
Michigan  48910 


Figure  1 :  Classifications  of  multiscale  methodologies 


properties.  Extrapolating  the  homogenized  response  to  non-linear  regions  or 
determining  composite  failure  will  rely  heavily  on  experimental  tests  to  characterize 
material  failure.  While  a  variety  of  failure  criteria  exist  at  the  homogenized  scale  of  the 
lamina,  or  even  the  full  laminate,  the  World  Wide  Failure  Exercises  (WWFE)  have 
shown  that  no  available  models  successfully  capture  UD  failure  under  complex  multi- 
axial  loading  states  or  for  multi-ply  laminates  [8].  Concurrent  models  utilize  domain 
decomposition,  whereby  detailed  microstructural  information  is  directly  coupled  at 
small  regions  of  interest  to  a  homogenized  macroscopic  domain  [9-11].  A  drawback 
of  these  methods  is  that  the  computational  cost  and  memory  requirements  grow 
significantly  with  increasing  size  of  the  region  of  interest  (i.e.  due  to  large  crack 
growth  or  diffuse  damage).  Additionally,  adaptive  remeshing  schemes  and 
transitionary  elements  are  required  to  account  for  unknown  a-priori  damage  paths  and 
the  coupling  of  multiple  lengths  scales  in  one  FE  mesh.  Semi-concurrent  models,  on 
the  other  hand,  preserve  microstructural  information  by  linking  each  integration  point 
at  the  macroscopic  level  to  a  unique  unit-cell  boundary  value  problem  (BVP)  at  the 
microscopic  level  [12-14],  The  BVP,  defined  on  a  representative  volume  element 
(RVE),  will  determine  the  macroscopic  response  of  the  macroscopic  integration  point. 
This  linking  between  the  two  scales  is  shown  schematically  in  Figure  1,  highlighting 
the  passing  of  information  between  the  weakly  coupled  scales.  Despite  the  increased 
computational  cost  of  semi-concurrent  schemes  over  sequential  methods,  the 
uncoupled  nature  of  each  BVP  allows  for  parallelization  of  the  necessary 
microstructural  solution. 

The  work  presented  utilizes  a  semi-concurrent,  nested  finite  element  (FE) 
methodology  similar  to  that  pioneered  by  Feyel  and  Chaboche  [15],  whereby  the  unit¬ 
cell  BVP  is  solved  using  the  FE  method.  In  contrast  to  other  semi-analytical  methods 
[16,17],  implementing  a  FE  solution  scheme  allows  for  the  implementation  of  explicit 
damage  modeling  through  cohesive  zone  modeling  [18,19]  or  through  the  extended- 
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Finite  Element  Method  (XFEM)  [3,20,21],  A  significant  limitation  to  current  semi- 
concurrent  models  is  the  available  boundary  conditions  with  which  one  can  drive  the 
unit-cell  BVP.  While  periodic  boundary  conditions  have  found  to  be  effective  over 
uniform  displacement  or  uniform  traction  boundary  conditions  for  determining  the 
effective  properties  of  an  RVE  [22,23],  recent  work  by  Coenen  et  al.  [24]  and  Aboudi 
and  Ryvkin  [25]  have  shown  that  semi-concurrent  models  with  periodic  boundary 
conditions  cannot  capture  the  true  material  behavior  in  the  presence  of  localized 
damage.  Although  a  variety  of  novel  boundary  conditions  have  been  presented  in 
literature  to  overcome  this  limitation  [24,26],  the  proposed  boundary  conditions  are 
modified  forms  of  periodicity  and  unable  to  account  for  strong  non-periodic  effects. 
The  proposed  work  focuses  on  developing  a  methodology  for  implementing  non¬ 
periodic  boundary  conditions  in  a  semi-concurrent  scheme  while  preserving  energy 
between  the  scales. 


MULTISCALE  FRAMEWORK 

The  two  main  components  in  the  semi-concurrent  multiscale  implementation 
are  the  localization  and  the  homogenization  rules.  Localization  refers  to  the  passing  of 
information  from  the  homogenized  global  scale  to  local  scale  of  the  RVE.  Conversely, 
homogenization  refers  to  the  determination  of  macroscopic  quantities  from  the  RVE, 
or  the  passing  of  information  from  RVE  to  the  global  integration  point.  In  a 
deformation  driven  FE  analysis,  the  localization  mles  provides  the  kinematic  coupling 
from  macroscopic  to  microscopic  domains.  Although  the  process  can  use  the 
macroscopic  strain  at  the  selected  integration  point,  the  work  that  follows  uses  the 
macroscopic  deformation  gradient,  P1* ,  as  shown  in  Figure  2.  The  macroscopic 
deformation  gradient  is  then  used  to  specify  necessary  boundary  conditions  for  the 
RVE’s  BVP.  Standard  procedures  assume  a  volume  average  relationship  between  the 
deformation  at  the  global  scale  and  that  of  the  RVE  shown  in  Equation  (1). 

F'*  =  V'laFI?dV'  (1) 

The  superscripts  M  and  m  represent  macroscopic  and  microscopic  quantities, 
respectively,  F  is  the  deformation  gradient  tensor,  and  I  b  is  the  reference  volume  of 
the  RVE.  The  way  in  which  the  macroscopic  deformation  gradient  provides  boundary 

conditions  (displacements,  or  M  )  on  the  RVE  are  provided  in  Equation  (2), 

u?  =  (Fff  -  l)ip  (2) 

where  X  is  the  vector  of  reference  configuration  coordinates.  The  prescription  of 
these  displacements  can  be  chosen  for  all  nodes  within  the  RVE,  all  nodes  along  the 
boundary,  or  restricted  to  the  vertices  in  the  case  of  periodicity.  The  homogenization 
procedure,  which  provides  the  macroscopic  stress  as  a  function  of  microscopic 
stresses,  also  involves  an  volume  averaging  relationship  and  is  provided  in  Equation 
(3). 

^=Ua^dVB  (3) 

The  stress  is  written  using  the  first  Piola-Kirkhoff  stress  tensor,  P  ,  for  convention,  as 
it  is  the  work  conjugate  to  the  deformation  gradient.  A  second  part  of  the 
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homogenization  process  is  determining  the  instantaneous  material  Jacobian,  J  ,  shown 
in  Figure  2. 
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Figure  2.  Overview  of  the  semi-concurrent  workflow  and  passing  of  information  between  the 
macroscale  (global)  and  microscale  (RVE/local). 


PROPOSED  Z-SCALAR  METHOD 

All  multiscale  schemes  are  required  to  preserve  the  Hill-Mandel  relation  which 
states  that  the  variation  of  work  at  the  global  scale  must  equal  the  volume  average  of 
the  variation  of  work  at  the  local  scale.  This  relation  is  shown  graphically  in  Figure  3 
as  well  as  in  Equation  (4), 

SWM  =  j'  SW™  dVB  (4) 

where  0  W  is  the  variation  of  work  at  a  particular  scale.  As  a  result  of  the  Hill-Mandel 
relation,  the  kinematic  and  constitutive  coupling  between  the  macroscopic  and 
microscopic  domain  are  constrained  to  satisfy  the  averaging  relations  provided  in 
equations  (1)  and  (3).  These  two  equations  state  the  volume  average  of  the  stress/strain 
field  variables  throughout  the  RVE  are  equal  to  the  corresponding  macroscopic 
variables.  Given  these  averaging  relations,  it  is  a  trivial  calculation  to  show  that  the 
Hill-Mandel  relation  is  satisfied,  and  these  calculations  are  given  for  various  types  of 
RVE  boundary  conditions  in  work  by  Kouznetsova  [12].  Coenen  et  al.  [24]  enforced 
boundary  conditions  satisfying  the  strain  averaging  relations,  then  proved  that  the 
volume  averaged  microscopic  stress  tensor  would  in  fact  satisfy  the  work  equivalence 
between  the  scales.  Equation  (5)  below  expresses  the  Hill-Mandel  relation  in  terms  of 
the  individual  deformation  gradient  and  stress  tensors. 
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Work  (Macro) 


Work  (Micro) 


Figure  3.  The  work,  or  more  precisely  the  variation  of  work,  at  both  the  macro  and  micro  scales  must 
be  maintained  for  the  multiscale  scheme  to  be  energetically  consistent. 

As  was  discussed  in  the  work  by  Coenen  et  al.  [24],  the  strain  averaging 
relation  (or  alternatively,  the  deformation  averaging  relation)  is  satisfied  when  the 
contributions  of  the  micro-fluctuation  on  the  volume  average  are  zero.  A  simplified 
form  of  this  statement  is  shown  in  Equation  (6), 


(6) 


which  is  derived  from  the  right  hand  side  of  Equation  (1)  and  using  Equation  (7).  It 
should  be  noted  that  the  strain  averaging  relation  is  being  written  with  respect  to  the 
deformation  gradient  rather  than  the  strain  tensor,  however  the  equations  still  hold. 
The  following  relation 


(7) 


describes  the  local,  current  coordinates  x  within  the  RVE  as  a  function  of  the 
macroscopic  deformation  gradient  and  a  function  of  the  microfluctuation  field,  w, 
which  is  dependent  on  the  microstructure.  This  superposition  of  a  macroscale 
influence  and  the  microfluctuation  field  is  graphically  represented  in  Figure  4. 


RVE  deformation 


Macroscale  influence 


Microscale  influence 
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Figure  4:  Schematic  of  the  RVE  deformation  as  a  superposition  of  macroscale  and  microscale 

influences. 


Figure  5 :  Illustration  of  the  localization  process,  which  has  deformed  the  RVE  to  produce  a  unique 
amount  of  elastic  strain  energy  defined  by  the  box  on  the  right. 


The  strain,  or  deformation,  averaging  relation  will  only  hold  if  last  integral 
term  in  (6)  is  equal  to  zero.  It  has  been  shown  in  literature  that  periodic  boundary 
conditions  satisfy  this  requirement.  Non-periodic,  non-uniform  boundary  conditions, 
however,  will  contain  a  microscale  influence  field  which  will  not  set  the  integral  term 
in  (6)  to  be  zero.  The  microscale  influence  field  will  be  an  unknown  solution, 
dependent  on  the  microstructure,  thus  the  form  of  the  macroscopic  stress  tensor  cannot 
be  solved  analytically  a-priori.  Instead,  it  is  assumed  in  this  work  that  the  volume 
average  of  the  RVE  stress  field  is  a  sufficient  approximation,  although  there  is  no 
guarantee  (due  to  the  nature  of  the  microfluctuation  field)  that  it  is  consistent  with 
Hill’s  energetic  conditions. 

In  the  first  iteration  of  developing  an  energetically  consistent  stress  tensor,  an 
equivalence  of  elastic  strain  energy  between  the  scales,  as  introduced  by  Hill  [6],  is 
directly  enforced.  First,  the  localization  process  is  performed  as  in  Figure  5  according 
to  a  given  macroscopic  quantity  of  deformation.  At  this  stage,  it  is  assumed  that  the 
localization  provided  a  sufficient  state  of  stress/strain  within  the  RVE  resulting  in  a 
given  amount  of  elastic  strain  energy,  shown  by  the  blue  box  in  Figure  5.  It  is  desired 
to  determine  the  macroscopic  stress  tensor  required  to  ensure  that  the  elastic  strain 
energy  at  the  global  level  matches  that  from  the  RVE.  The  formulation  to  follow 
assumes  that  the  microconstituents  are  linear  elastic  with  cohesive  based  interactions 
between  the  fibers  and  matrix.  Due  to  the  presence  of  elastic  softening  from  the 
failing/failed  cohesive  surfaces,  the  elastic  strain  energy  at  the  macroscopic  level  will 
be  computed  according  to  the  diagram  in  Figure  6-(b),  rather  than  the  purely  elastic 
case  shown  in  Figure  6-(a).  The  2D  plane  strain  problem,  to  which  this  work  is 
currently  focused,  requires  that  the  elastic  strain  energy  is  computed  from  the 
contributions  of  stress  and  strain  in  the  3  degrees  of  freedom  in  the  problem  (11,  12, 
and  22  directions),  shown  in  Figure  7.  The  condition  of  equivalent  elastic  strain  energy 
between  the  scales  enforces  that  the  total  energy  shown  in  Figure  7  should  be  equal  to 
that  found  during  the  localization  process  in  Figure  5.  Thus,  this  equality  provides  a 
constraint  with  which  we  can  determine  the  form  of  the  macroscopic  stress.  Although 
there  are  an  infinite  number  of  stress  tensors,  given  known  values  of  deformation 
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(strain)  at  the  macroscopic  level,  that  will  satisfy  the  energy  equivalence,  it  is  hereby 
assumed  that  the  macroscopic  stress  state  can  be  approximated  as  being  proportional 
to  the  volume  average  of  the  RVE  stresses. 


Figure  6:  Stress-strain  diagrams  highlighting  the  elastic  strain  energy  for  two  non-linear  material 
responses  where  (a)  the  material  is  purely  elastic  and  (b)  the  material  has  non-linearity  due  to  elastic 

softening. 


This  assumption  of  proportionality  introduces  a  scaling  parameter,  z,  which 
will  be  found  by  enforcing  energy  equivalence  between  the  global  and  local  scale.  The 
computation  of  the  macroscopic  stress  tensor  will  now  be  found  using 


ff-T  =  z  * 


<=zxl. 


G[?dV0 


(8) 


where  &*M  is  the  macroscopic  Cauchy  stress  tensor  computed  from  volume  averaging 
the  microscale  Cauchy  stress  tensors,  cr''1 ,  in  the  RVE.  The  energy  balance  with 
respect  to  the  elastic  strain  energy  at  the  macro  and  micro  scales  can  be  written 


— 

** s  train 


+  j«)= z .  gsJK  + 


_L _ cM 

T-  2  t23 


where  left  hand  side  represents  the  elastic  strain  energy  at  the  microlevel,  ^strain , 
from  the  RVE  localization  shown  in  Figure  5,  while  the  right  hand  side  represents  the 
elastic  strain  energy  computed  using  the  volume  averaged  stresses.  The  z  scalar 
parameter  is  computed  using 


z  = 


strain 


f\.Af  _*Af  1l 

Vi  1  ^1 1  +elSffl!  +  £^7H  ) 


(10) 


where  it  is  a  function  of  the  microscopic  elastic  strain  energy,  macroscopic  strains,  and 
volume  averaged  RVE  stresses.  This  scalar  parameter  represents  this  energy  based 
constitutive  coupling  between  the  microscopic  and  macroscopic  scales.  The  next 
subsection  will  discuss  the  numerical  implementation  of  this  novel  concept  into  the 
semi-concurrent  scheme. 
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Figure  7 :  The  total  amount  of  elastic  strain  energy  at  the  macroscopic  level  for  the  2D  plane  strain 
problem  is  a  summation  of  the  elastic  strain  energy  resulting  from  the  three  degrees  of  freedom. 


NUMERICAL  IMPLEMENTATION 

The  semi-concurrent  scheme  as  well  as  the  proposed  z-scalar  methodology 
was  implemented  numerically  utilizing  Python  scripting  to  invoke  the  nested  FE 
solution  within  the  commercial  FE  software  ABAQUS.  To  reduce  initial  overhead  in 
implementing  the  multiscale  framework  the  current  model  was  restricted  to  2D,  plane 
strain  analysis.  The  iterative  algorithm  is  presented  below  in  Figure  8.  At  each 
incremental  step  in  the  analysis,  the  user  material  defined  subroutine  (UMAT)  was 
utilized  to  perform  the  communication  between  Python  and  the  ABAQUS  solver.  The 
left-hand  side  of  Figure  8  highlights  the  localization  process  which  involves  passing  of 
the  macroscopic  deformation  gradient  from  the  UMAT  to  the  custom  Python  script 
which  then  modifies  the  boundary  conditions  to  a  unit-cell,  or  RVE,  ABAQUS  input 
file.  The  unit-cell  job  is  submitted  and  once  completed  is  post-processed  according  to 
the  right-hand  side  of  Figure  8.  As  a  test  case  for  implementing  non-periodic  boundary 
conditions,  a  choice  was  made  to  implement  boundary  conditions  equivalent  to  that 
shown  in  Figure  9-(a).  These  test  boundary  conditions  contains  a  free  edge  at  one 
boundary,  corresponding  to  a  macro  scale  with  similar  boundaries.  The  symmetric 
edge  at  the  left  of  the  RVE  was  an  approximation  chosen  arbitrarily. 
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LOCALIZATION 


UMAT: 

l  i  mat  subroutine  outputs  the 
deformation  gradient  F*'  to  <t  ,t*t  tile 


Python: 

R.?<iils  Fr-‘ 

Opens  Unit  Cell  input  file 
Modifies  tlie  hoimddty  conditions 
of  tlie  unit-cell  input  file 

Submits  the  job  in  ABACUS 


HOMOGENIZATION 


UMAT: 

Redds  in  die  mdcroscopic  stress 
Redds  in  the  cut  rent  Jdcohidit 
Submits  to  ARAOI  IS 


Python: 

Redds  unit  cell  .COB  file 
Computes**  from  volume  averaging 
C dlcu Idles  tlie  anient  Jdcobidn P 


Figure  8:  Workflow  of  the  semi-concurrent  computational  homogenization  scheme. 


Figure  9:  Schematic  of  (a)  the  test  boundary  conditions,  which  include  the  free  edge,  and  (b)  the 
standard  periodic  boundary  conditions.  Note:  That  P,  F,  S  in  the  above  diagram  represent  periodic, 

free  and  symmetric  edge,  respectively. 


The  implementation  of  the  free-edge  boundary  conditions  results  in  a  different 
set  of  localization  rules  than  those  used  for  periodic  boundary  conditions  given  as 

u?  =  F*fX?  (11) 

The  localization  in  equation  (11)  is  identical  to  that  in  (2)  except  that  the  microscale 
displacements  are  only  prescribed  on  the  vertices,  v,  of  the  RVE.  Elsewhere,  non¬ 
vertex  nodes  are  constrained  to  obey  periodicity  according  to  the  formulation  of  Van 
der  Sluis  et  al.  [22]  and  using  *Equation  keywords  in  ABAQUS.  In  the  case  of  the 
boundary  conditions  shown  in  Figure  9-(a),  the  localization  shown  in  Equation  (11)  is 
prescribed  for  vertices  1  and  4.  For  vertices  2  and  3,  displacements  are  only  prescribed 
in  the  2-direction  to  account  for  the  free-edge  requirement  which  must  remain  traction 
free  in  the  1-direction.  Additionally,  all  nodes  lying  on  the  symmetric  edge  (S*)  are 
prescribed  displacements  in  the  1 -direction  and  left  undefined  in  the  2-direction. 
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1.  Reads  the  EVOL,  or  current  element  volume  for  all  elements 

2.  Extracts  the  element  stresses  within  the  unit-cell  RVE 

3.  Computes  the  volume  average  of  the  stresses  based  on  the  EVOL 
values 

4.  Using  the  volume  averaged  stresses  and  Equation  (10),  z  is  calculated 

5.  The  macroscopic  stress  tensor  is  exported  to  an  external  “.csv”  file 


Figure  10.  Listing  of  the  database  post-processing  steps  performed  in  the  custom  Python  script. 


When  implementing  periodic  boundary  conditions,  step  4  in  Figure  10  is  omitted. 
Aside  from  the  macroscopic  stress  tensor,  the  material  Jacobian  needs  to  be 
determined  for  the  subsequent  increment  in  the  global  FE  analysis  and  is  a  required 
output  of  the  UMAT  within  ABAQUS  Standard.  The  construction  of  the  Jacobian 
from  the  perturbation  steps  is  summarized  in  Figure  1 1 .  For  a  given  perturbation  step, 
a  component  of  strain  is  set  to  unity,  while  all  other  are  set  to  zero.  Thus,  the  resulting 
stresses  computed  as  a  result  provide  a  given  column  in  the  jacobian  stiffness  matrix. 
The  stresses  are  computed  using  the  same  procedure  as  was  done  for  the  macroscopic 
stress  tensor,  utilizing  the  scaling  parameter  when  necessary  (e.g.  when  non-periodic 
localization  is  employed).  The  J33  component  must  be  found  externally  using  standard 
micromechanics,  or  basic  homogenization  techniques.  Assuming  the  current  type  of 
RVE,  little  damage  will  be  accumulating  in  the  fiber  direction,  thus  this  component 
should  remain  a  constant  through  the  entire  2D  analysis.  The  remaining  components  in 
the  Jx3  column  are  found  assuming  symmetry  of  the  Jacobian.  Once  the  Jacobian  is 
computed  and  exported  to  a  “.csv”  file,  the  Python  script  is  completed.  The  UMAT 
will  then  read  in  the  macroscopic  stress  tensor  and  Jacobian  which  outputted  back  into 
the  ABAQUS  simulation  at  the  global  scale.  When  periodic  boundary  conditions  are 
employed,  the  perturbation  steps  are  executed  immediately  after  the  stress  analysis 
within  the  same  ABAQUS  job.  For  the  case  of  the  test  boundary  conditions  in  Figure 
9-(a),  the  operation  must  be  performed  as  a  separate  unit-cell  job  after  the  completion 
of  the  stress  analysis.  Additionally,  the  boundary  conditions  are  modified  to  remove 
the  free-edge  and  are  shown  in  Figure  12.  The  right-hand  side  of  the  RVE  is  perturbed 
using  equation  (11)  for  all  nodes  in  both  the  1-  and  2-directions.  The  process  for 
determining  the  Jacobian  during  the  post-processing  steps  is  listed  in  Figure  13.  An 
overall  summary  of  the  discussed  workflow  with  the  Fortran  subroutine  and  Python 
script  are  highlighted  in  Figure  14. 
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Set  =1 


Figure  11.  The  Jacobian  matrix,  J,  is  shown  with  the  necessary  uni-strain  components  for 
determining  a  particular  column.  The  33  component  is  determined  externally. 


Figure  12.  Boundary  conditions  used  for  the  perturbation  steps  during  non-periodic  analysis. 


6.  Extract  the  nodal  coordinates  at  the  end  of  the  localization  step  for  all 
perturbed  nodes 

7.  Use  Equation  (14)  to  determine  BC’s  for  the  three  perturbation  steps 
Note:  These  steps  are  defined  as  linear  perturbation  steps  in  ABAQUS 

8.  Generate  the  necessary  input  file  job  and  submit  (Perturb.inp) 

9.  Compute  the  differential  macroscopic  stress  tensor  associated  with  each 
perturbation 

10.  Build  the  Jacobian  and  export  the  resultant  to  an  external  “.csv”  file 


Figure  13.  Listing  of  the  post-processing  steps 


Verification 

The  first  step  in  verifying  the  modifications  coded  into  the  semi-concurrent 
scheme  was  to  check  the  validity  and  robustness  of  the  Jacobian  updating  procedure 
using  the  boundary  conditions  discussed  in  the  previous  section.  This  test  was  carried 
out  by  comparing  computed  elastic  moduli  (using  Mathematica  and  the  exported 
Jacobian  matrices  from  Python)  to  those  predicted  from  purely  periodic  boundary 
conditions.  A  number  of  test  cases  were  examined  including  a  homogeneous  RVE,  a 
single  fiber  RVE,  and  multi-fiber  RVE.  For  the  fiber/matrix  RVE’s  the  fiber  volume 
fraction  was  45%.  Table  I  presents  a  comparison  of  the  reference  periodic  values  with 
those  computed  using  the  boundary  conditions  shown  in  Figure  12.  For  all  three  test 
cases,  the  non-periodic  perturbations  provided  reasonable  estimations  of  the  material 
constants.  It  should  be  noted,  however,  that  due  to  the  uniform  displacement  boundary 
conditions  on  the  right-hand  side  of  the  RVE  caused  the  elastic  moduli  to  be  slighlty 
over-predictive. 
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TABLE  I.  Computed  elastic  constants  for  periodic  and  non-periodic  perturbation  steps. 


Material 

Constant 

Homogeneous  Solid 

Single  Fiber  RVE 
(Em  =3.5  GPa, 

Ef  =  22GPa) 

25-fiber  RVE 
(Em  =  3.5  GPa, 

Ef  =  22GPa) 

Periodic 

Non- 

Periodic 

Periodic 

Non- 

Periodic 

Periodic 

Non- 

Periodic 

Ex  (GPa) 

100.00 

100.000 

20.897 

21.257 

20.247 

20.455 

Ey  (GPa) 

100.00 

100.000 

20.897 

21.137 

20.378 

20.419 

vxy 

0.300 

0.300 

0.294 

0.295 

0.330641 

0.328911 

Gxy  (GPa) 

38.462 

38.462 

6.438 

6.901 

7.207075 

7.363027 
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Figure  14.  Complete  workflow  of  the  numerical  implementation  of  the  non-periodic,  semi- 

concurrent  scheme. 
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The  next  step  in  the  verification  process  was  to  confirm  that  the  z-scalar 
method  preserved  energy  between  the  global  and  local  scales  even  when  using  non¬ 
periodic  boundary  conditions.  A  single,  reduced-order  quadrilateral  element  (CPE4R) 
at  the  global  scale  was  placed  under  uniaxial  tensile  loading.  The  element’s  integration 
point  was  coupled  to  a  25  multi-fiber  RVE  with  fibers  placed  randomly  within  the 
boundary  of  the  RVE.  Both  periodic  as  well  as  the  non-periodic  free  edge  boundary 
conditions  were  employed  in  two  separate  iterations  and  the  global  elements  were 
stretched  to  1%  strain.  As  was  assumed  in  the  formulation  of  the  z-scalar  method,  all 
constituents  are  linear-elastic,  however,  a  cohesive  interface  following  a  bilinear 
traction  separation  exists  between  the  fiber  and  matrix.  The  necessary  parameter  for 
the  cohesive  interface  are  the  cohesive  strength,  &cahe  ,  and  the  critical  displacement, 
dcrit.  The  cohesive  strength  defines  the  peak  traction  value  in  the  traction-separation 
law,  and  the  critical  displacement  specifies  the  traction-free  separation  distance.  The 
cohesive  parameters  are  chosen  to  be  representative  of  a  relatively  weak  interface.  The 
constituent  and  interface  properties  are  provided  in  Table  II. 


TABLE  II.  Constituent  properties  for  the  multi-fiber  RVE  verification  study. 


E  (GPa) 

V 

^cohe 

dcrit 

(MPa) 

(m) 

Fiber 

22 

0.3 

- 

- 

Matrix 

3.5 

0.4 

- 

- 

Interface 

10E7 

- 

50 

0.005 

It  was  found  during  the  verification  studies  that  the  ETOTAL,  or  total  energy 
history  variable,  of  the  system  was  non-zero,  a  phenomenon  that  was  indicative  of 
incorrect  energy  computations.  Figure  15  shows  all  the  energy  history  variables 
outputted  by  ABAQUS.  Further  investigations  revealed  that  the  negative  ETOTAL, 
highlighted  by  the  dashed  box  in  Figure  15,  were  a  result  of  the  contributions  of  the 
cohesive  surface  to  the  internal  strain  energy  of  the  system.  When  using  cohesive 
interactions  rather  than  cohesive  elements,  the  strain  energy  present  in  the  separation 
of  the  fiber/matrix  interface  was  not  getting  added  to  the  total  strain  energy  of  the 
RVE  system.  Although  the  contribution  of  the  cohesive  interface  to  the  strain  energy 
of  the  system  can  be  reduced  by  increasing  the  cohesive  stiffness,  even  a  sufficiently 
large  value  of  10E7  GPa  had  a  measurable  contribution  as  shown  in  Figure  15. 
Equation  (12)  shows  the  energy  balance  equations  of  relevant  energy  sources  in  the 
simulations. 
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Due  to  the  neglecting  of  cohesive  strain  energy,  the  Estra iai ,  and  consequently  the 
internal,  energy  terms  were  being  under  reported  resulting  in  a  negative  ETOTAL 
result. 
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Figure  15.  Energy  history  variable  outputs  for  the  16  fiber  RVE  under  periodic  boundary  conditions 

subject  to  uniaxial  tensile  loading. 
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■^“Periodic 
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Figure  16.  Percent  difference  in  external  work  between  the  macro  and  micro  scale  work  for  the  16  fiber 
RVE  for  periodic  and  non-periodic  boundary  conditions. 


As  a  consequence  of  the  negative  ETOTAL  values,  the  computation  of  the  z- 
scalar  parameter  in  Equation  (10)  from  the  microstructural  strain  energy  had  to  be 
modified  as  in 


(13) 


pm  —  p  _  P 

"strata  "strata  ^ total 

The  values  E strain  and  E total  in  equation  (13)  were  extracted  from  ABAQUS  from 
the  ALLE  and  ETOTAL  history  variables. 
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The  results  of  the  z-scalar  verification  are  presented  in  Figure  16,  plotting  the 
percent  difference  in  external  work  computed  by  ABAQUS  between  the  global  and 
local  scales  for  the  two  boundary  conditions  cases  (periodic  and  non-periodic  with  a 
free  edge).  Both  the  global  and  local  external  work  energies  were  found  to  be  nearly 
identical  for  the  two  test  cases,  as  seen  from  the  relatively  small  values  (<0.4%)  in 
Figure  16.  Coincidentally,  the  two  test  cases  had  similar  external  work  values 
regardless  of  the  differences  in  applied  boundary  conditions  on  the  level  of  the  RVE. 
For  the  case  of  the  non-periodic  boundary  conditions,  the  z-scalar  was  computed  to  be 
a  value  slightly  less  than  unity.  For  example,  at  the  last  increment  of  the  macroscale 
analysis  (-1%  strain)  the  z-scalar  parameter  was  computed  as  0.984.  It  was 
hypothesized  that  this  small  variation  from  unity  was  a  function  of  the  test  boundary 
conditions  and  a  lack  of  significant  differences  in  the  RVE  response  for  the  periodic 
and  non-periodic  BC’s.  Since  periodic  boundary  conditions  have  been  proven  to  be 
energetically  consistent,  the  close  agreement  between  the  macroscale  periodic  and 
non-periodic  external  work  in  Figure  16  would  suggest  that  the  implemented  z-scalar 
parameter  ensures  similar  consistency.  Also,  the  small  variations  seen  even  with  the 
periodic  boundary  conditions  were  hypothesized  to  be  numerical  construct  within 
ABAQUS,  rather  than  an  imbalance  of  work  between  the  two  scales. 

APPLICATION 

A  case  study  was  undertaken  to  show  an  application  of  z-scalar  energy  based 
method.  The  2D  RVE  used  for  the  energy  verifications  and  discussed  in  the  previous 
sections  was  used  in  the  following  case  study.  The  objective  was  to  investigate  the 
effects  of  a  free  edge  boundary  at  the  macroscopic  level  on  the  constitutive  response 
provided  by  the  RVE.  The  boundary  conditions  for  the  RVE  were  the  test,  free-edge 
boundary  conditions  proposed  in  Figure  9-(a).  The  RVE’s  constituent  properties  were 
those  reported  in  Table  II.  Again,  damage  was  restricted  to  fiber/matrix  debonding. 
Fibers  were  placed  into  the  RVE  using  an  iterative  random  placement  script  written  in 
Python,  for  which  the  highest  achievable  fiber  volume  fraction  was  45%.  Macroscopic 
loading  was  uniaxial  tension  as  in  the  diagram  on  Figure  16.  Differences  between 
periodic  and  non-periodic  RVE  boundary  conditions  were  evaluated  according  to  the 
macroscopic  response  (evaluated  by  plotting  the  elastic  strain  energy  vs  strain)  as  well 
as  the  elastic  moduli  (Ex,  Ey).  Previous  research  works  [27,28]  have  shown  that  RVE 
size  (number  of  fibers)  affects  the  macroscopic  response,  and  the  appropriate  RVE 
size  is  dependent  on  the  constituent  properties  of  the  RVE.  Other  works  have  plotted 
the  variance  in  macroscopic  (stress-strain)  response  over  a  variety  of  boundary 
conditions,  none  of  which  included  a  free-edge  [23,24,26].  In  this  case  study,  the  RVE 
size  is  varied  from  4  fibers  to  36  fibers  with  multiple  iterations  (3-4)  of  random  fiber 
placement. 

Figure  17  plots  the  elastic  strain  energy  versus  macroscopic  strain  for  the  4 
fiber  RVE,  across  several  random  placement  iterations,  as  well  as  the  comparisons 
between  periodic  and  non-periodic  (free)  boundary  conditions.  The  curves  are  labeled 
according  to  XPerY  or  XFreeY,  where  X  is  the  RVE  size  and  Y  is  the  iteration 
number.  Dashed  curves  represent  the  periodic  simulations,  and  solid  curves  represent 
those  with  the  non-periodic  boundary.  The  sudden  plateaus  seen  at  -0.6%  strain,  or 
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dip  in  the  case  of  iterations  2  &  3,  are  a  result  of  the  cohesive  failure.  The  curves  were 
as  expected  from  trend  seen  in  the  elastic  energy  curve  shown  in  Figure  15. 
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Figure  17.  Specific  elastic  strain  energy  vs  macroscopic  strain  for  the  4-fiber  RVE.  Dashed  lines 
represent  the  periodic  results,  solid  lines  those  of  the  free  edge  simulations,  and  the  color  coding 

corresponds  to  a  given  RVE  iteration. 
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Figure  18.  Macroscopic  (2-direction)  stress  vs  strain  for  the  4-fiber  RVE. 

Figures  18-20  plot  the  macroscopic  stress  in  the  loading  (22)  direction  versus 
macroscopic  strain.  These  plots  provide  the  RVE’s  macroscopic  material  response 
under  the  uniaxial  tensile  loading.  For  the  4  RVE  cases  in  Figure  18,  only  one  case 
(iteration  3)  had  a  significantly  different  response  between  the  two  BC’s.  A  similar 
observation  can  be  made  from  Figure  17.  Figure  18  also  highlighted  the  significant 
variation  in  material  behavior  (cohesive  failure)  between  the  4  different  RVE 
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iterations  as  a  function  of  the  random  placement.  Iterations  1  &  4  had  significantly 
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Figure  19.  Macroscopic  (2-direction)  stress  vs  strain  for  the  16-fiber  RVE. 
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Figure  20.  Macroscopic  (22  direction)  stress  vs  strain  for  the  36-fiber  RVE. 

lower  peak  stress  values  (pre-cohesive  failure)  and  a  more  drastic  drop  in  stress  before 
reloading.  Figure  19  shows  the  results  from  the  16  fiber  test.  The  simulation  for  the 
16Per3  terminated  early  due  to  convergence  issues  present  after  cohesive  failure  and  is 
indicated  on  the  figure.  The  first  16  fiber  RVE  iteration  had  a  nearly  identical 
macroscopic  response  between  periodic  and  non-periodic,  while  the  second  iteration 
had  significant  softer  post-peak  response  in  the  periodic  BC  case.  For  both  of  these 
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iterations,  the  periodic  response  was  “softer”,  indicative  of  increased  cohesive  damage 
as  a  result  of  the  BC.  The  significant  differences  in  iterations  1  and  2  between  the 
periodic  and  non-periodic  counterparts  could  be  attributed  to  the  RVE  fiber  placement 
and  effect  of  the  free-edge.  The  RVE  for  iteration  1,  shown  at  the  top  left  of  Figure  19, 
had  fewer  fibers  in  the  proximity  of  the  right-hand  edge  versus  the  RVE  for  iteration 
2,  shown  at  the  bottom  left  of  the  figure.  The  periodic  BC  caused  cohesive  failure  and 
separation  at  the  fiber  highlighted  by  the  black  arrow  on  Figure  19,  whereas  the  free 
edge  did  not  cause  failure/separation.  The  36  fiber  RVE  results  presented  in  Figure  20 
highlighted  two  main  points.  First,  as  expected  from  literature,  the  variation  in 
material  macroscopic  response  decreased  across  the  random  fiber  iterations.  Second, 
the  periodic  boundary  conditions  suffered  from  numerical  convergence  issues 
resulting  from  premature  failure  in  the  cohesive  surfaces,  highlighted  on  the  figure. 
This  happened  for  two  of  the  three  periodic  cases;  however,  the  premature  failure  did 
not  affect  the  material’s  tensile  strength  as  was  with  the  premature  failure  of  the 
16Per3  case.  The  free-edge  simulations,  however,  were  able  to  reach  the  final 
macroscopic  loading  of  1%  strain  with  no  issues.  Iteration  1,  which  completed  to  1% 
for  both  BC’s,  showed  little  differences  in  the  macroscopic  response.  The  similarity  of 
periodic  and  non-periodic  response  for  iteration  1  can  again  be  attributed  to  the 
proximity  of  fibers  to  the  right  edge  of  the  RVE,  shown  at  the  top-right  of  Figure  19. 
Additionally,  extracting  the  initial  modulus  (slope)  from  the  across  the  periodic 
(dashed)  and  non-periodic  (solid)  curves  in  Figures  18-20.  First,  there  was  a  decrease 
in  variation  of  the  initial  modulus  within  increasing  RVE  size.  The  initial  modulus 
converged  for  the  36  fiber  case  (6.7  GPa),  and  was  slightly  softer  than  that  predicted 
for  the  4  and  16  fiber  cases  (7. 1-6.8  GPa). 


CONCLUSION 

The  developed  energy  based  coupling  using  the  z-scalar  method  allows  one  to 
investigate  the  effect  of  non-periodic  RVE  boundary  conditions  at  the  microscale  on 
the  global  response  of  the  composite.  The  introduction  of  a  z-scalar  parameter  allows 
for  the  preservation  of  elastic  strain  energy  between  the  macroscopic  and  microscopic 
domain  and  is  independent  of  the  type  of  boundary  conditions  employed  on  the  RVE. 
The  current  implementation,  however,  is  limited  to  a  2D  RVE  domain,  and  the 
formulation  of  the  z-scalar  method  restricts  the  RVE  to  one  for  which  elastic  softening 
is  the  cause  of  macroscopic  nonlinearity  (in  this  case,  cohesive  failure). 

Despite  these  current  limitations,  the  robustness  of  the  proposed  methodology 
was  shown  through  a  practical  investigation  of  the  effect  of  macroscopic  free-edges  on 
the  macroscopic  response  of  a  2D  fiber/matrix  RVE  with  cohesive  interfaces.  The 
results  presented  suggest  that  the  overall  macroscopic  response  was  largely  affected  by 
fiber  proximity  to  the  free-edge  at  the  microscopic  level.  A  larger  number  of  fibers 
close  to  the  free-edge  result  in  greater  damage  evolution  from  periodic  BC’s  versus 
non-periodic,  free  edge  BC’s.  The  results  suggest  that  higher  fiber  volume  fraction 
RVE’s,  creating  an  abundance  of  fibers  near  the  free-edge,  would  increase  the  effect 
of  the  free-edge  to  the  damage  evolution.  Namely,  a  stiffer  post-peak  response  would 
result  from  accounting  for  a  macroscopic  free-edge. 

Future  work  will  develop  the  z-scalar  methodology  into  a  work-based 
equivalence,  rather  than  elastic  strain  energy,  to  allow  for  plastic  deformation  and/or 
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crack  growth  within  the  RVE.  Applying  a  work-based  form  of  the  z-scalar  method  and 
novel  BC’s  will  allow  for  the  exploration  RVE  crack  initiation  and  modeling  using 
XFEM.  Appropriate  localization  schemes  will  be  developed  to  allow  for  deformations 
conducive  for  crack  opening. 
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